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ABSTRACT 

In 3D image/video acquisition, different views are often cap¬ 
tured with varying noise levels across the views. In this paper, 
we propose a graph-based image enhancement technique that 
uses a higher quality view to enhance a degraded view. A 
depth map is utilized as auxiliary information to match the 
perspectives of the two views. Our method performs graph- 
based filtering of the noisy image by directly computing a 
projection of the image to be filtered onto a lower dimen¬ 
sional Krylov subspace of the graph Laplacian. We discuss 
two graph spectral denoising methods: first using Chebyshev 
polynomials, and second using iterations of the conjugate 
gradient algorithm. Our framework generalizes previously 
known polynomial graph filters, and we demonstrate through 
numerical simulations that our proposed technique produces 
subjectively cleaner images with about 1-3 dB improvement 
in PSNR over existing polynomial graph filters. 

Index Terms — 3D images, image denoising, graph- 
based filtering, Krylov, Chebyshev, conjugate gradient 

1. INTRODUCTION 

In 3D video applications, multiple views are often presented 
at different quality levels due to various acquisition and com¬ 
pression conditions applied to the source signals. For ex¬ 
ample, changes of brightness or color may be produced by 
imaging sensors and circuitry of stereo cameras or even from 
shot noise. Traditional image enhancement techniques may 
be used to improve the quality of the degraded views. More¬ 
over, depth map compression has become an integral part of 
emerging 3D video formats, e.g., 3D-AVC and 3D-HEVC. 
Hence, it is desirable to exploit depth information to enhance 
a low quality view. 

Graph signal processing tools have recently been applied 
to classical image processing tasks EEiiiiiiia. For exam¬ 
ple, a typical interpolation problem was studied using spectral 
graph theory in O, where the upsampling problem was for¬ 
mulated as a regularized least squares problem. Later, this 
approach was extended to depth image upsampling in (41, 
and slight benefits of graph spectral domain processing were 


demonstrated. Recently, a similar graph based method to en¬ 
hance noisy stereo images was investigated in 13 utilizing 
the depth information to generate the guide image. However, 
these methods suffered from high complexity due to their re¬ 
quirement of computing a full eigenvalue decomposition on 
very high dimensional data. 

An alternative approach that avoids full eigen- 
decompositions was adopted in O [6| [71 where the graph 
spectral filter was approximated by degree-/c polynomials. 
The filtering operation was then performed by applying the 
same polynomial as a function of the graph Laplacian matrix 
in the pixel domain. 

We follow the same line of thought and define a general 
framework for graph spectral filtering by directly computing 
the projection of the desired filtered signal onto an order- 
{k 1) Krylov subspace for the graph Laplacian. We dis¬ 
cuss how image denoising can be realized by low-pass fil¬ 
tering with a degree-/c Chebyshev polynomial with a prede¬ 
fined stop band. We then propose a parameter free filtering 
approach that utilizes k iterations of the conjugate gradient 
(CG) algorithm to achieve the same denoising objective. We 
demonstrate through numerical experiments that our filters 
better suppress noise while maintaining image details com¬ 
pared to previous methods. 

The remainder of the paper is organized as follows. Sec¬ 
tion 2 introduces the notation, gives a basic background on 
graph-based image processing, and reviews prior-art methods. 
Section 3 constructs a simple 4-connected graph using a guide 
image warped by 3D projection, and then proposes Cheby¬ 
shev and CG filters to handle the image denoising problem. 
Section 4 describes our numerical experiments and analyses 
the results, demonstrating the improved performance of our 
filters. Section 5 concludes the work. 

2. GRAPH BASED IMAGE PROCESSING 
2.1. Basics of images on graphs 

In graph signal processing m, an undirected graph G = 
(L, E) consists of a collection of vertices, also called nodes, 
V = {1,2, ...,A} connected by a set of edges E = 
{{i, j,Wij)},i, j G V, where {i,j,Wij) denotes an edge be- 


tween nodes i and j associated with a weight Wij >0. A 
degree di of a node i is the sum of edge weights connected 
to the node i. For image processing applications, a pixel may 
be treated as a node in a graph, while the edge weight may be 
viewed as a measure of similarity of the pixels. 

An adjacency matrix W of the graph is a symmetric 
N X N matrix having entries Wij > 0, and a diagonal degree 
matrix is D := diagjdi, ^ 2 ,... A graph Laplacian ma¬ 

trix L := D — W is a positive semi-definite matrix, thus ad¬ 
mitting an eigendecomposition L = UAU^, where U is an 
orthogonal matrix with columns forming an orthonormal set 
of eigenvectors, and A = diag{ Ai,..., Aat} is a matrix made 
of corresponding eigenvalues, all real. The smallest eigen¬ 
value of the matrix L is always zero. 

Alternatively, a normalized symmetric Laplacian matrix 
is defined as Ls := D-V2LD-1/2, and is a positive semi- 
definite matrix. Hence, it admits an eigendecomposition 
Ls = UsA^Us^, where Us is an orthogonal matrix con¬ 
structed from an orthonormal set of eigenvectors and As is 
its corresponding diagonal eigenvalue matrix. The smallest 
eigenvalue of the matrix Ls is always zero, while the largest 
eigenvalue is conveniently bounded by two. 

The eigenvalues and eigenvectors of the Laplacian ma¬ 
trixes provide a spectral interpretation of graph signals, where 
the eigenvalues can be treated as graph Fourier frequencies, 
and the eigenvectors as generalized Fourier modes. 

A graph design can be associated to an underlying con¬ 
ventional image filter, traditionally used as one of the filter 
testing benchmarks. Specifically, for an input image the 
conventionally filtered output image can be written as, 

Xont = D“^Wx^n = Xin “ D“^Lx^n- (1) 

The original L or the symmetric normalized Laplacian Ls 
can serve as a foundation for image processing in the graph 
spectral domain. The peculiarity of the symmetric normal¬ 
ized Laplacian Ls is that it has all its eigenvalues located on 
the interval [0, 2] and that it requires pre- and post-processing 
of the images, i.e., the signal in vertex domain needs to be 
normalized x^^ = D 2 before applying the filter, and de- 
normalized back Xont = D“ 2 Xont after being filtered. 

Remark: In the rest of the paper, we assume, for certainty, 
that the symmetric normalized Laplacian Ls is always used, 
and we drop the subscript “S”, denoting the symmetric nor¬ 
malized Laplacian as L = UAU^ and calling it simply the 
“Laplacian” hereafter. For example, conventional filter Q in 
the new notation turns into Xout = ^in — Lx^^. 

Graph Spectral Filtering (GSF) FL can be designed for im¬ 
age processing purposes in the graph spectral domain, where 
FL is 3 . diagonal matrix, typically given as = h{A), where 
/i(A) is a real valued function of a real variable A, determin¬ 
ing the filter. The corresponding graph filter H in the vertex 


(a) 


(b) 


Fig. 1: Graph structure (a) used in 111; (b) proposed in Section 
S.llof this paper. 


domain can be expressed as, 

H = /i(L) = UHU^. (2) 

For example, taking /i(A) = A, we obtain the graph Laplacian 
h{L) = L itself that constitutes a high pass filtering operator 
attenuating low frequency coefficients and boosting high fre¬ 
quency coefficients of a signal. 

In E and m the filtered signal is computed by solving 
the following regularized least squares problem, 

X* = arg min h\x-hf + ^||/lRx||^ (3) 

2 2 

where b := x^^ is the noisy signal. Hr is a penalty high pass 
graph filter, and p is a regularization parameter. We refer to 
this formulation as graph based joint bilateral filter (GBJBF). 
The above problem admits an analytic solution of the form, 
X* = U(I + pnn^)-^V^h = (I + pH|)-ib. 

Choosing FLn = A results in the penalty filter Hr = L. 
Consequently, the GBJBF filter E is given by, 

^gbjbf(A) = (1 + pX ^)~^. (4) 

In E, the filter function /igbjbf(*) 0 is approximated by a 
truncated degree-/c Chebyshev polynomial expansion, which 
we denote by /ik-POLY(')' 

We call a filter “polynomial,” if h{X) is a polynomial func¬ 
tion in A. For example, GBJBF filter /igbjbf(-) ^ is not 
polynomial, while its approximation /ik-POLY(-) is. The con¬ 
ventional filter Q is given by h{X) = 1 — A, which is a first 
degree polynomial of A, i.e., is also polynomial. We note that 
application of the conventional filter via Q does not require 
knowledge of the eigendecomposition of the Laplacian L, in¬ 
volving instead only a calculation of the product of the matrix 
L and a given vector. 

This is a common feature of polynomial filters, including 
the filter given by /ik-POLY(') E . making polynomial filtering 
computationally attractive. If the underlying polynomial is 
given via its roots or by a recursive formula, its action on an 
image can be implemented, using only the calculation of the 
product of the matrix L and a vector. 

2.2. Graph-based image denoising 

A graph structure needs to be defined before setting and ap¬ 
plying graph based approaches for image denoising. In one 





























example, a 4-connected simple graph structure is as shown in 
Fig. [2 (a). Each pixel is connected only to its four immediate 
neighbors in the pixel space. All pixels within the image or 
an image slice share a single graph and would be processed 
within one graph spectral domain. In Section [Tllwe propose 
an adapted version of the graph structure in Fig.[^(b) to han¬ 
dle 3D image denoising with help from depth information. 

After the connection structure in the graph is defined, a 
weight needs to be assigned for each graph edge. A well- 
known approach is to assign bilateral weights such as in m 
and m, where the weights Wij are defined by, 

m, = (5) 

The first exponential term is a spatial distance penalty (pi 
refers to the pixel’s spatial location) and the second exponen¬ 
tial term is an intensity distance penalty (xin refers to the in¬ 
tensity value from a guide image). We note that for the 4- 
connected simple graph structure as shown in Fig. (a) the 
first exponential term is a constant. The guide image is the 
image itself in ID and a warped image in Q. Once the joint 
bilateral graph is constructed, the underlying filter as defined 
in 0 is referenced as Joint Bilateral Filter (JBF) hereinafter. 

3. PROPOSED CHEBYSHEV AND CG GRAPH 
FILTERING 

3.1. Proposed graph structure 

We use a graph structure similar to the 4-connected graph 
shown in Fig. ^h). The edge weights of the links are de¬ 
fined as in <0 with guided image being the warped image. 
A standard depth image based rendering (DIBR) process is 
employed to generate the guidance image as described below. 

For each pixel location i = [u,v] in the current view, a 
corresponding location i' = [u', v'] can be determined based 
on the camera parameters and the pin hole camera model. 
Hence, Xin{i) in the warped image is obtained from x[^{i') 
in the other high quality view. Note that although [u^ v] is 
always at a pixel at integer location, [u' ^v'] may point to a 
subpixel location. As a result, subpixel interpolation in the 
high quality view needs be performed to maintain accuracy. 
In this paper, we use the 8- or 7- tap interpolation filter de¬ 
fined in H.265/HEVC video coding standard (21. 

Occlusions may occur near foreground objects during the 
warping process, which are marked as hole areas. A typical 
hole filling algorithm would propagate the background pix¬ 
els into the holes using inpainting algorithms. Since the filled 
holes are often unreliable, we do not use the filled pixels as 
guidance to construct the graph. Instead, we propose to sim¬ 
ply avoid any links to or from a hole pixel. As shown in 
Fig.[T](b), as there is a hole pixel below pixel pi, there are 


three links left for this pixel. More complex graph connec¬ 
tions are subject to future research. 

A conventional filter corresponding to the underlying fil¬ 
ter of the above graph is herein referenced as JBF, since joint 
bilateral weights are used. 

3.2. Graph filtering via subspace projection 

The graph filtering problem can be viewed in the context of a 
general task of applying a function /i(L) of the graph Lapla- 
cian L to an input signal thereafter for brevity denoted by 
b, such that, 

X* = /i(L)b. (6) 

For image denoising, the goal is to suppress high frequency 
noise in which case the graph filter /i(L) is a low pass filter. 

For general functions /i(L), e.g., such as the heat kernel 
it becomes intractable to exactly compute the action of 
/i(L) on b via the eigendecomposition of L if the dimension¬ 
ality of b is large. To overcome this difficulty, we propose 
instead to setup and evaluate the projection of x* onto an ap¬ 
propriate low dimensional subspace JC. Denoting by P^; the 
projector onto the subspace /C, the projection x^ of the solu¬ 
tion X* of 0 in the subspace 1C is given by, 

x^ = Pjch{L)h. (7) 

If the subspace 1C is L-invariant, i.e., L/C C 1C then 
= /i(P/cLPx:), and thus 'P)ch{'L)h can be effi¬ 
ciently computed for an arbitrary function /i('), e.g., utilizing 
a basis of the subspace 1C made of eigenvectors of L. 

The optimal choice of the L-invariant subspace 1C is to 
apply a principal component analysis to the matrix h{Jj) and 
select the subspace 1C spanned by the principal eigenvectors, 
corresponding to the largest by absolute value eigenvalues of 
This choice implies that 'P]ch{'L) ~ /i(L) is the best 
small-rank approximation. Since the matrices h(L) and L 
share the same eigenvectors, one can actually compute eigen¬ 
vectors of L corresponding to the largest by absolute value 
eigenvalues of /i(L). If the shape of the function h{') is 
known in advance, one may be able to select the regions of 
interest in the spectrum of L for eigenvector computations. 
We note that this construction of the subspace 1C does not 
explicitly depend on x^^ = b. It may offer a dramatic com¬ 
putational cost improvement, compared to the full eigenvalue 
decomposition of /i(L) or L, if an efficient iterative method 
is used to compute the basis of the subspace 1C. 

In this paper, we are interested specifically in low-pass 
filters, so we concentrate on using the Krylov subspace pro¬ 
jection technique li uni, i.e., we approximate the solution x* 
by a vector in the order-{k + 1) Krylov subspace, 

1C = span{b, Lb,..., L^b}. (8) 

The Krylov subspace is known to approximate well eigen¬ 
vectors corresponding to extreme eigenvalues, thus it should 




be an appropriate choice for high- and low-pass filters. Al¬ 
though Krylov subspace I®. is not normally exactly L- 
invariant, we can still attempt using hiVjcLVjc) as the only 
practically available, within the constraints of polynomial fil¬ 
tering, replacement for V]ch{L). Moreover, Krylov subspace 
@ is built starting with the initial image = b, therefore it 
only takes into account the L-spectral modes actually present 
in the image. 

An implementation of a filter using the projection on the 
degree k Krylov subspace can be especially simple if the fil¬ 
ter function h{') itself is a polynomial of degree k or smaller, 
h{X) = Pk{X), since = P/cP/c(L)b = p/e(L)b, where the 
equality holds since the subspace 1C contains all possible vec¬ 
tors p/c(L)b for any polynomial Pk{') of degree k or smaller, 
i.e. Pk (L)b G 1C. Computational gain is achieved when the 
matrix L is sparse and the degree k is small. 

3.3. Chebyshev polynomial graph spectral denoising 

In the graph-based image denoising problem, we wish to ap¬ 
ply a low pass filter to the graph spectrum of the noisy im¬ 
age in order to remove high frequency noise. Restricting our¬ 
selves to polynomial filtering in this work, e.g., due to compu¬ 
tational considerations, we can setup a problem of designing 
optimal polynomial low pass filters. Our first optimal polyno¬ 
mial low pass filter is based on Chebyshev polynomials. 

Specifically, we propose to use a degree k Chebyshev 
polynomial /ik-CHEB defined over the interval [0, 2] with a stop 
band extending from I G (0, 2) to 2. Since we define the graph 
spectrum in terms of the eigenspace of the symmetric normal¬ 
ized Laplacian L, all the eigenvalues of L lie in the interval 
[0, 2]. The construction of a Chebyshev polynomial is easily 
obtained by computing the roots of the degree k Chebyshev 
polynomial r{i) = cos {7r{2i — l)/2k) for i = 1... k. over 
the interval [—1,1], then shifting the roots to the interval [/, 2] 
via linear transformation to obtain the roots of the polyno¬ 
mial /ik-CHEB, and scaling the polynomial using vq such that 
^k-CHEB(O) = 1. This results in formula. 



Fig. 2: Spectral responses of the graph filters, (a) Differ¬ 
ent graph filters: JBF, GBJBF, k-POLY O, and proposed k- 
CG; (b) Proposed CG filter with different sequences: Kendo, 
Pozan_Street, and Undo_Dancer. 



Fig. 3: Spectral response with at different iterations, (a) k- 
POLY lO; (b) k-CG in this paper. 


graph spectral filter by choosing h(L) = L^, i.e., setting 
X* = L^b. However, in contrast to the approach of m , where 
^gbjbe(') is approximated by a truncated Chebyshev polyno¬ 
mial expansion, we do not attempt to approximate the inverse 
function by a polynomial. Instead, we use the subspace pro¬ 
jection technique, formulating the graph filtering problem as 
a constrained quadratic program, 

Xp = argmin x^Lx — 2x^f, 

_ xg£ (9) 

/C = x^ + span {f — Lx^,..., L^(f — Lx^)} , 


^k-i 


k-CHEB I 


i=l 


(A) = pk(A) = n (^ ) = ^0 n ; 




and we can compute x^ by evaluating x* = riyC~^ — Lx*“^ 
iteratively for i = 1,..., /c, where x^ = rob. 

Chebyshev polynomials are minimax optimal, uniformly 
suppressing all spectral components on the interval [/, 2] and 
growing faster than any other polynomial outside of [/, 2]. The 
stop band frequency I remains a design parameter that needs 
to be set prior to filtering. In the next subsection, we propose 
a variational parameter-free method. 


3.4. Conjugate Gradient method for Krylov denoising 

Since L is a high-pass filter, its Moore-Penrose pseudoinverse 
Lt is a low-pass filter, so we propose designing a low-pass 


where the initial approximation x^ and the vector f remain to 
be chosen. The solution to problem (|^ can be computed by 
running k iterations of the CG method ifTTIl . We denote by 
hk-cG the k-step CG filter with x^ = f = b, in which case 
iC C 1C —the order-(/c+1) Krylov subspace defined in We 
denote by /ik-CGO the k-stop CG filter with x^ = b and f = 0. 
Yet another choice x^ = 0 and f = b results in IC equal to 
the order-k Krylov subspace span{b, Lb,..., L^“^b}. 

Fig. (a) compares the spectral response h{X) of three 
iterations of the proposed 3-CG approach with the benchmark 
methods: JBF, GBJBF, and 3-POLY—a degree-3 polynomial 
^ 3 -POLY approximation of GBJBF as in im. The CG algorithm 
adaptively adjusts the spectral response based on the input 
image as shown in Fig.[^(b). Fig.[^(a) and (b) plot /ik-POLY 
and hk-cG and illustrate the variation in spectral response over 
different iterations for k-POLY and k-CG with /c = 1, 2, 3. 

































Table 1: Denoised image PSNR results (dB) 


Seqs 

Kendo 

Poznan.Street 

Undo-Dancer 

JBF 

32.53 

32.78 

31.90 

3-POLY 

32.28 

33.05 

31.97 

3-CG 

35.76 

34.49 

31.91 

3-CHEB 

35.67 

34.35 

31.92 




(a) Original Image. (b) Noisy Image. 



(c) Filtered Image, 3-POLY. (d) Filtered Image, 3-CG. 
Fig. 4: Comparison of filtered images. 


4. EXPERIMENTS AND DISCUSSIONS 

We evaluate the performance of the proposed approaches for 
graph filtering on stereo plus depth 3D video denoising. In 
our experimental setup, one view is captured at high qual¬ 
ity while the other view is degraded with Gaussian noise. 
The depth map associated with the stereo video is avail¬ 
able. Three videos with varying characteristics are chosen: 
Kendo @ 1024 x 768 and Poznan_Street are captured from 
natural scenes with estimated depth maps; Undo_Dancer @ 
1920 X 1080 is a synthetic sequence generated from 3D mod¬ 
els with ground truth depth maps ca. To filter pixels in the 
hole area, we simply apply a 3 x 3 median filter. 


4.1. Performance evaluations 


We compare the denoising performance of four approaches 


using the same graph structure as proposed in Section 3.1 


The first benchmark method is JBF, which shares the same 
weights as the underlying graph. The second benchmark 


Fig. 5: Comparison of spectral response of the proposed 3- 
CHEB and 3-CG filters. 


method is 3-POLY, which is a degree-3 polynomial approx¬ 
imation to the regularization solution 0 of GBJBF with pa¬ 
rameter p = 2. We then run our degree-3 Chebyshev filter 
(3-CHEB) and 3 iterations of our CG approach (3-CG) and 
compare the denoised image to the benchmark techniques. 

The denoised image PSNR results are summarized in Ta¬ 
ble [T] On the one hand, all methods demonstrate similar de¬ 
noising performance for the Undo_Dancer sequence, where 
the PSNR differences vary between 0.01 - 0.07 dB. This re¬ 
sult might be explained by the lack of complex texture in this 
synthetic sequence. 

On the other hand, we observe that our CG approach re¬ 
sults in a significant improvement of 3.23 dB and 1.44 dB 
in PSNR over the best benchmark techniques for sequences 
Kendo and Poznan_Street, respectively. The quality improve¬ 
ment is also verified in subjective evaluations. Fig.|^(c) and 
(d) are the filtered images from 3-POLY and 3-CG. Notice 
that the CG approach preserves texture details better than 3- 
POLY, e.g., plate numbers resulted from the CG approach are 
easier to recognize. The content adaptation feature of CG is 
illustrated by various spectral response curves for different 
sequences in Fig.|^(b). 

The Chebyshev filter that we design has a stop band at I = 
0.5. Notice that the denoised image PSNR is almost identical 
to the CG PSNR. That can be explained by our choice of /, 
which is close to the effective stop band of our CG filter k- 
CGO. We plot examples of the spectral response of our k-CG, 
k-CGO, and k-CHEB filters in Fig.j^ 

4.2. Patch-based processing 

In order to reduce memory requirements in implementing the 
above filtering schemes, we use 64 x 64 disjoint patches of 
the image that are filtered independently. Moreover, the patch 
based approach is favorable for parallel implementation and 
could be synchronized with the coding unit structures used in 
a typical video coding framework such as H.265/HEVC. 

A subjective evaluation of the patch based processing ap¬ 
proach shows that GBJBF sometimes exhibits block artifacts 
across patch boundaries. Fig.[^(c) illustrates these block ar¬ 
tifacts as witnessed by the difference between the noisy and 
denoised images. On the other hand, our CG approach does 
not suffer from blocky denoising as shown in Fig.[^(d). The 
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(a) Original Image. (b) Mask Map. 



(c) Diff. Image, 3-POLY. (d) Diff. Image, 3-CG. 

Fig. 6: Block artifacts along patch boundaries. 


only block viewed in the difference image of the CG approach 
is due to the switch between the median filter applied to the 
hole pixels, shown in Fig.|^(b), and the CG technique. 

Finally, the polynomial degree k, which we set equal to 3 
for 64 X 64 patches, tends to increase with larger patch sizes. 
We note that for k-POLY, the degree needs to be determined 
beforehand so as to determine the polynomial coefficients and 
cannot be incrementally increased without recalculating the 
coefficients and restarting the filtering process. 

In contrast, the k-CHEB approach allows incrementally 
increasing the degree k, although only in such a way that the 
higher degree polynomials share the same roots as the previ¬ 
ously computed lower degree polynomials, which is possible 
for some degrees which are powers of 2 and 3, or their prod¬ 
ucts. Moreover, the k-CG technique allows for a seamless in¬ 
cremental increase in the degree of the polynomial by simply 
increasing the number of iterations. 


5. CONCLUSIONS AND FUTURE WORK 

We propose Chebyshev and conjugate gradient filters for 
graph based image denoising. Our approach performs graph- 
based filtering on the noisy image by directly computing the 
projection of the desired filtered image onto a lower dimen¬ 
sional Krylov subspace for the graph Laplacian using novel 
Chebyshev or conjugate gradient polynomial graph filters. 
We demonstrate through numerical simulations that our pro¬ 
posed technique produces subjectively cleaner images with 
about 1-3 dB improvement in PSNR over existing polynomial 
graph filters. The proposed approaches are not limited to 3D 
image denoising and could be used for improving other types 
of images, provided that guidance images are constructed un¬ 
der proper contexts, which is subject to future work. 
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